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Abstract 

Clustering methods have been used extensively in computer vision and pattern 
recognition. Fuzzy clustering has been shown to be advantageous over crisp (or traditional) 
clustering in that total commitment of a vector to a given class is not required at each 
iteration. Recently fuzzy clustering methods have shown spectacular ability to detect not 
only hypervolume clusters, but also clusters which are actually "thin shells", i.e., curves 
and surfaces. Most analytic fuzzy clustering approaches are derived from Bezdek's Fuzzy 
C-Means (FCM) algorithm. The FCM uses the probabilistic constraint that the 
memberships of a data point across classes sum to one. This constraint was used to 
generate the membership update equations for an iterative algorithm. Unfortunately, the 
memberships resulting from FCM and its derivatives do not correspond to the intuitive 
concept of degree of belonging, and moreover, the algorithms have considerable trouble in 
noisy environments. Recently, we cast the clustering problem into the framework of 
possibility theory. Our approach was radically different from the existing clustering 
methods in that the resulting partition of the data can be interpreted as a possibilistic 
partition, and the membership values may be interpreted as degrees ot possibility of the 
points belonging to the classes. We constructed an appropriate objective function whose 



minimum will characterize a good possibilistic partition of the data, and we derived the 
membership and prototype update equations from necessary conditions for minimization of 
our criterion function. In this paper, we show the ability of this approach to detect linear 
and quartic curves in the presence of considerable noise. 


1 


Research performed for NASA/JSC through a subcontract from the RICIS Center at the 


University of Houston - Clear Lake 


I. Introduction 


Clustering has long been a popular approach to unsupervised pattern recognition. It 
has become more attractive with the connection to neural networks, and with the increased 
attention to fuzzy clustering. In fact, recent advances in fuzzy clustering have shown 
spectacular ability to detect not only hypervolume clusters, but also clusters which are 
actually "thin shells", i.e., curves and surfaces [1-7]. One of the major factors that 
influences the determination of appropriate groups of points is the "distance measure" 
chosen for the problem at hand. Fuzzy clustering has been shown to be advantageous over 
crisp (or traditional) clustering in that total commitment of a vector to a given class is not 
required at each iteration. 

Boundary detection and surface approximation are important components of 
intermediate-level vision. They are the first step in solving problems such as object 
recognition and orientation estimation. Recently, it has been shown that these problems can 
be viewed as clustering problems with appropriate distance measures and prototypes [1-7]. 
Dave's Fuzzy C Shells (FCS) algorithm [2] and the Fuzzy Adaptive C-Shells (FACS) 
algorithm [7] have proven to be successful in detecting clusters that can be described by 
circular arcs, or more generally by elliptical shapes. Unfortunately, these algorithms are 
computationally rather intensive since they involve the solution of coupled nonlinear 
equations for the shell (prototype) parameters. These algorithms also assume that the 
number of clusters are known. To overcome these drawbacks we recently proposed a 
computationally simpler Fuzzy C Spherical Shells (FCSS) algorithm [6] for clustering 
hyperspherical shells and suggested an efficient algorithm to determine the number of 
clusters when this is not known. We also proposed the Fuzzy C Quadric Shells (FCQS) 
algorithm [5] which can detect more general quadric shapes. One problem with the FCQS 
algorithm is that it uses the algebraic distance, which is highly nonlinear. This results in 
unsatisfactory performance when the data is not very "clean" [7], Finally, none of the 


algorithms can handle situations in which the clusters include lines/planes and there is much 
noise. In [frl, we addressed those issues in a new approach called Piano-Quadric 
Clustering. In this paper, we show how that algorithm, coupled with our new possibilistic 
clustering, can accurately find linear and quadric curves in the presence of noise. 

Most analytic fuzzy clustering approaches are derived from Bezdek's Fuzzy C- 
Means (FCM) algorithm [9]. The FCM uses the probabilistic constraint that the 
memberships of a data point across classes must sum to one. This constraint came from 
generalizing a crisp C-Partition of a data set, and was used to generate the membership 
update equations for an iterative algorithm. These equations emerge as necessary conditions 
for a global minimum of a least-squares type of criterion function. Unfortunately, the 
resulting memberships do not represent one's intuitive notion of degrees of belonging, i. 
e„ they do not represent degrees of "typicality" or "possibility". 

There is another important motivation for using possibilistic memberships. Like all 
unsupervised techniques, clustering (crisp or fuzzy) suffers from the presence of noise in 
the data. Since most distance functions are geometric in nature, noise points, which are 
often quite distant from the primary clusters, can drastically influence the estimates of the 
class prototypes, and hence, the final clustering. Fuzzy methods ameliorate this problem 
when the number of classes is greater than one, since the noise points tend to have 
somewhat smaller membership values in all the classes. However, this difficulty still 
remains in the fuzzy case, since the memberships of unrepresentative (or noise) points can 
still be significantly high. In fact, if there is only one real cluster present in the data, there is 
essentially no difference between the crisp and fuzzy methods. 

On the other hand, if a set of feature vectors is thought of as the domain of 
discourse for a collection of independent fuzzy subsets, then there should be no constraint 
on the sum of the memberships. The only real constraint is that the assignments do really 
represent fuzzy membership values, i.e., they must lie in the interval [0,1]. In [10], we cast 



the clustering problem into the framework of possibility theory. We briefly review this 
approach, and show it's superiority to recognize shapes from noisy and incomplete data. 

II. Possibilistic Clustering Algorithms 

The original FCM formulation minimizes the objective function given by 
C N C 

J C L,U ) = Z £ (M ii ) m d\- , subject to . Z u.. = 1 for ally. ( 1 ) 

In (1), L = (Aj,...,A c ) is a C-tuple of prototypes, dy is the distance of feature point Xj to 

cluster A,-, N is the total number of feature vectors, C is the number of classes, and U = 
iMijl is a C xN matrix called the fuzzy C-partition matrix [9] satisfying the following 

conditions: 
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Here, fly is the grade of membership of the feature point Xj in cluster A z , and m e [ 1 ,°°) 
is a weighting exponent called the fuzzifier. In what follows, A f - will also be used to denote 
the ith cluster, since it contains all of the parameters that define the prototype of the cluster. 

Simply relaxing the constraint in (1) produces the trivial solution, i. e., the criterion 
function is minimized by assigning all memberships to zero. Clearly, one would like the 
memberships for representative feature points to be as high as possible, while 
unrepresentative points should have low membership in all clusters. This is an approach 
consistent with possibility theory [11]. The objective function which satisfies our 
requirements may be formulated as: 
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where 77 / are suitable positive numbers. The first term demands that the distances from the 

feature vectors to the prototypes be as low as possible, whereas the second term forces the 
11- to be as large as possible, thus avoiding the trivial solution. The following theorem, 

proved in [9], gives necessary conditions for minimization, hence, providing the basis for 
an iterative algorithm. 

Theorem: 

Suppose that X = {xj, ..., x^,} is a set of feature vectors, L = (Aj,...,A c ) is a 
C-tuple of prototypes, ctj is the distance of feature point Xj to the cluster prototype A ( -, (/ 
= 1, ..., C; j = 1, ..., N), and U = [/i~] is a C xN matrix of possibilistic membership 

values. Then U may be a global minimum for J m (L,U) only if = 

The necessary conditions on the prototypes are identical to the corresponding conditions in 
the FCM and its derivatives. 

Thus, in each iteration, the updated value of fly depends only on the distance of x,-. 
from A,-, which is an intuitively pleasing result. The membership of a point in a cluster 
should be determined solely by how far it is from the prototype of the class, and should not 
be coupled to its location with respect to other classes. The updating of the prototypes 
depends on the distance measure chosen, and will proceed exactly the same way as in the 
case of the FCM algorithm and its derivatives. 

The value of r\ , determines the distance at which the membership value of a point in 
a cluster becomes 0.5 (i. e., "the 3 dB point"). Thus, it needs to be chosen depending on 
the desired "bandwidth" of the possibility (membership) distribution for each cluster. This 



value could be the same for all clusters, if all clusters are expected to be similar. In general, 
it is desirable that 77, relates to the overall size and shape of cluster A,. Also, it is to be 
noted that 77, determines the relative degree to which the second term in the objective 

function is important compared to the first. If the two terras are to be weighted roughly 
equally, then 77, should be of the order of . In practice we find that the following 


definition works best. 

N 



( 3 ) 


This choice makes 77, the average fuzzy intra-cluster distance of cluster A,. The value of 77, 
can be fixed for all iterations, or it may be varied in each iteration. When 77, is varied in 
each iteration, care must be exercised, since it may lead to instabilities. Our experience 
shows that the final clustering is quite insensitive to large (an order of magnitude) 
variations in the values of 77,. 


III. The Possibilistic C Piano-Quadric Shells Algorithm 

Suppose that we are given a second degree curve Aj characterized by a prototype 

vector 

T 

Pi = \Pil'Pi 2 Pir\ 

to which it is desired to fit points Xj obtained through the application of some edge 
detection algorithm. If a point x has coordinates (x j, ... , x n ]. then let 
2 2 2 

q = [.tj, x 2 x n ,X[X 2 >*(/!- l)*/j’*l’ x 2 x n> U 7 • 


When the exact (geometric) distance has no closed-form solution, one of the methods 
suggested in the literature is to use what is known as the "approximate distance" which is 


the first-order approximation of the exact distance. It is easy to show [12] that the 
approximate distance of a point from a curve is given by 
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where ^^Qjj is the gradient of the distance functional 
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evaluated at xj . In (4) the matrix Dj is simply the Jacobian of q evaluated at xj . 

One can easily reformulate the quadric shell clustering algorithm with as the 

underlying distance measure. It was shown in [8] that the solution to the parameter 
estimation problem is given by the generalized eigenvector problem 

Pi Pi = li Gi pi , (6) 


where 


N 

Fi = . 



Gi = . Dj Dj T , 

which can be converted to the standard eigenvector problem if the matrix Gi is not rank- 
deficient. Unfortunately this is not the case. In fact, the last row of Dj is always [0, . . , 
,0]. Equation (6) can still be solved using other techniques that use the modified Cholesky 
decomposition [13], and the solution is computationally quite inexpensive when the feature 
space is 2-D or 3-D. Another advantage of this constraint is that it can also fit lines and 


planes in addition to quadrics. Our experimental results show that the resulting algorithm, 
which we call the Possibilistic C Piano-Quadric Shells (PCPQS) algorithm, is quite robust 
in the presence of poorly defined boundaries (i. e., when the edge points are somewhat 
scattered around the ideal boundary curve in the 2-D case and when the range values are not 
very accurate in the 3-D case). It is also very immune to impulse noise and outliers. Of 
course, if the type of curves required are restricted to a single type, e.g., lines, or circles, 
or ellipses, simpler algorithms can be used with possibilistic updates, as will be seen. 

IV. Determination of Number of Clusters 

The number of clusters C is not known a priori in some pattern recognition 
applications and most computer vision applications. When the number of clusters is 
unknown, one method to determine this number is to perform clustering for a range of C 
values, and pick the C value for which a suitable validity measure is minimized (or 
maximized) [14]. However this method is rather tedious, especially when the number of 
clusters is large. Also, in our experiments, we found that the C value obtained this way 
may not be optimum. This is because when C is large, the clustering algorithm sometimes 
converges to a local minimum of the objective function, and this may result in a bad value 
for the validity of the clustering, even though the value of C is correct. Moreover, when C 
is greater than the optimum number, the algorithm may split a single shell cluster into more 
than one cluster, and yet achieve a good value for the overall validity. To overcome these 
problems, we proposed in [8] an alternative Unsupervised C Shell Clustering algorithm 
which is computationally more efficient, since it does not perform the clustering for an 
entire range of C values. 

Our proposed method progressively clusters the data starting with an overspecified 
number Cmax of clusters. Initially, the FCPQS algorithm is run with C=Cmax ■ After the 
algorithm converges, spurious clusters (with low validity) are eliminated; compatible 


clusters are merged; and points assigned to clusters with good validity are temporarily 
removed from the data set to reduce computations. The FCPQS algorithm is invoked again 
with the remaining feature points. The above procedure is repeated until no more 
elimination, merging, or removing occurs, or until C= 1. 

V. Examples of Possibilistic Clustering for Shape Recognition 

Figures 1 and 2 show the detection of a circular "fractal edge" from a 
synthetically generated image. Figure 1(a) is the original composite fractal image; figure 
1(b) shows what a gray-scale edge operator finds (or doesn't find); figure 1(c) is the output 
of the horizontal fractal edge operator; with figure 1(d) giving the maximum overall 
response of the fractal operators in four directions. Figure 2(a) depicts the (noisy) 
thresholded and thinned result from figure 1(d). Figure 2(b) gives the final prototype found 
by the FPQCS (which, since there is only one cluster present, is the same as the crisp 
version). Note how the presence of noise distorts the final prototype. Figure 2(c) shows 
the possibilistic algorithm output, which is superimposed on the original image in figure 
2(d). The results of the PPQCS algorithm are virtually unaffected by noise. Several 
examples comparing crisp, fuzzy and possibilistic versions of clustering can be found in 
[ 6 , 8 , 10 ], 

Figure 3 depicts the algorithm applied to the image of a model of the Space Shuttle. 
Figure 3(a) is the original image. Figure 3(b) gives the output of a typical edge operator. 
Note that, due to the rather poor quality of the original image, the edges found both noisy 
and incomplete. This data was then input into the possibilistic piano-quadric clustering 
algorithm. Figure 3(c) gives the eight complete prototypes which were found after running 
the algorithm. Finally, figure 39(d) displays the prototype drawn only where sufficient 
edges points exist. 

VI. Conclusions 



In this paper, we demonstrated how our new possibilistic approach to objective- 
function-based clustering coupled with our piano - quadric shells algorithm can recognize 
first and second degree shapes from incomplete and noisy edge data. This approach is 
superior to both crisp and fuzzy clustering, as well as to traditional methods such as the 
Hough Transform. Extensions of this approach to other classes of shapes is currently 
underway. 
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Figure 1. Detection of a fractal circular edge. 

(a) Upper Left. Original fractal composite image. 

(b) Upper Right. Output of gray scale edge operator. 

(c) Lower Left. Output of "horizontal" fractal edge operator. 

(d) Lower Right. Results of Maximum magnitude c: outputs of four directions of fractal operators. 






Figure 2. Recognition of circular boundary. 

(a) Upper Left. Figure 1(d) thresholded and thinned. 

(b) Upper Right. Circular prototype found by fuzzy (or crisp) clustering. 

(c) Lower Left. Circular prototype found by possibilisiic clustering. 

(d) Lower Right. Possibilistic prototype superimposed on original image. 



Figure 3. Recognition of Shuttle model boundaries. 

(a) Upper Left. Original Shuttle image. 

(b) Upper Right. Incomplete and noisy edges found by edge operator. 

(c) Lower Left. Prototypes found by Possibilistic Piano-Quadric clustering. 

(d) Lower Right. Possibilistic prototypes superimposed drawn where there is sufficient edge 

information. 








Pose Estimation Using Possibilistic Clustering 

In the Third Quarter report, we described how the Unsupervised C Quadric Shells 
(UCQS) algorithm could be used to estimate the pose of the shuttle. The shuttle's image is 
taken from the back so that the exhaust nozzles and the back edges of the three wings are 
apparent. Given an original unrotated image, the exhaust nozzles can be parametrized by 
three circles, and the three wings can be parametrized by three straight lines. These 
parameters are easily determined by the UCQS algorithm. As the shuttle rotates, the shape 
of the nozzles will change from circles to ellipses, so will the orientation of the straight 
lines representing the three wings. The UCQS algorithm is used in order to cluster this 
edge image and determine the parameters of the ellipses and lines. Finally, these parameters 
can be used to solve for the translation and rotation parameters, as long as the translation is 
made in the image plane. In fact, depth information can also be derived from the change in 
the size of the nozzles. 

We also consider the case where only line information is available. Once again, our 
new possibilistic piano-quadric clustering approach is used to detect and recognize the 
linear segments. In what follows, derivation of pose parameters is given for both the case 
where three corresponding line segments have been identified, and where one circle and 
one line have been matched. 



POSE ESTIMATION: 


The 3-D object attitude in space can be determined from a single perspective image. 
Dhome et al [1] developed a method to solve for the three-dimensional attitude of an object 
based on the perspective projection of three image lines. Krishnapuram & Casasent [2] 
developed a method for determining two of the three rotation angles necessary to describe 
an object attitude in 3-D space from a single perspective projection of one circle. 

I. Determination of The Attitude of One Object From Three Lines: 

The perspective projection of a point Pi = (Xi, Yi, Zi) on an image is the point pi 
= (xi, yi, zi) = (Xi f/Z, Yj f/Z, f). Let li be an image line characterized by a vector vi = 
(aj, bj, 0) and a point pi = (xj, yi, f). li is the perspective projection of a space line Lj. 
Therefore it lies in the "interpretation plane" containing the origin of the coordinate system 
O and the image line li. The normal Nj to this plane is perpendicular to vj and the vector 
Opi. Thus Nj = vi ¥ Opi = (bj f, -aj f, di) T , where di = aj yj - bi xj is the Euclidean 
distance between the center of the image and line li. If Vi = (Ai, Bi, Ci) T is the director 
vector of the space line Lj, then it must be orthogonal to Nj, hence Vi . Ni = 0 implying 
that: 

(Aj, Bi, Ci) T . (bi, -aj, dj/f) = 0 (1) 

Consider three object lines in 3-D space Loi, i = 1, .... 3 defined in a model reference frame 
(Som)- The director vector of Loi is Voi = (Aoi, Boi. Q)i) T . When the object is rotated in 3- 
D space, the lines Loi are rotated into lines L 31 . Therefore 

(A 3 i, B 3 i, C3i) T = Rapy ( A 0i> Boi, Q)i) T ( 2 ) 

where RaPy is the rotation matrix. 

The perspective projections of lines L3j are the lines lQi- Equation (1) becomes 


( 3 ) 


(A3i, B 3 i, C3i) T . (boi, -a()i, d()i/f) T = 

Rapy (Aoi, Boi, Coi) T . (boi, -aoi. d()i/f) T = 0 

where i = 1,..., 3 and a, p, and y are the unknown rotation angles about x, y, and z axes 
respectively. Solving this system of equations is too complicated. A specially defined 
model coordinate system (Sim) and a corresponding viewer coordinate system (Siv) can 
be used to simplify the problem [1], With these coordinate systems, only two rotation 
angles a and P need to be determined, i.e. the system of equations (3) can be reduced to 
two equations and two unknowns. First, a is found by iteratively solving an 8^ order 
equation. Then P is solved for by substitution. When the three lines are coplanar, or when 
they form a junction, the 8^ order equation reduces to a 4th order equation. 

II Determination of the Attitude of an Object From a Circle and a Line: 


Given a circular curve on the x-y plane, and an x' y’ view of this curve in a 
different coordinate system x' y' z'. The two frames (x, y, z) and (x - . y\ z') are related by 
a homogenous transformation T, such that 
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A circle of radius r on the xy plane is described by : 

( o 2 2 

lx" + y = r 

lz = 0 


(4) 

(5) 


In the (x\ y'. z') frame, equations (4) & (5) become 



(tl 1 x' + t!2y'+tl3z') 2 + (t 2 lx' + t 22 y’+t 23 z') 2 = r 2 


( 6 ) 

(7) 


t3lx’ + t32y’+t33z’ = 0 


Substituting z' in terms of x' and y' from equation (7) into equation (6) yields the equation 
for the 2-D projection of the 3-D circular curve onto an arbitrary x' y' plane. Making use of 
the fact that the columns of T are mutually orthogonal unit vectors, we obtain 


(1 + 
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( 8 ) 


This is the equation of an ellipse in the (x\ y', z') frame. If the parameters of this ellipse 
are known, equation (8) can be solved for the transformation parameters t3i, t32, and t33. 
The transformation matrix T can be written as a function of the rotation angles a, p, and y. 


T = 


cosy cosP 

siny cosP 

-sinP 

0 


cosy sinP sina - siny cosa cosy sinP cosa + siny sina 0 
siny sinP sina + cosy cosa siny sinP cosa - cosy sina 0 

cosP sina cosP cosa 0 

0 0 1 


Having already solved for t3i, t32, and t33 , a and b can be easily determined from the 3 rd 
row of T. 

In order to determine the 3 rd angle y, a line can be used in addition to the circle. In this case 
the two rotation angles a and P can be determined as discussed previously. Knowing these 
two angles, equation (3) with i = 1 (since we have only one line) becomes simple to solve, 
since the only unknown is y. 
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